!========================================================
     Program main
     implicit none

!-----------------------------------------------
!     Variable declarations
!-----------------------------------------------
      integer nfiles, mtime, nmonth, nregion, mT
      parameter (mT=70)

      character*100 casename 
      character*100 cfile, cfile1
      character*50 ctemp

!-----------------------------------------------------
!     other variables
!-----------------------------------------------------
      integer ii, ij, ik, j, irec, ikk, irec_tot, itime, imon, ifile
      integer iregion, idx1, idx2, iix, iT
      integer ncount_rh(2, 100, mT)
      real*4 dist_rh(2, 100, mT), rh(100)
      real*4 total_ncount(2, mT)
      real*4 rhi1, rhi2, temp1, temp2
      real*4 a, b, c, tc
      integer idx_25, idx_50, idx_75
      real*4  ni_mean
      real*4  accum_pdf(100)
      real*4  tk, nice, rei


      ncount_rh(:, :, : ) = 0.0
      dist_rh(:, :, :) =  0.0
      do idx1=1, 100 
!  1.0e0 to 1.0e3 
        rh(idx1) = 1.0e0 * 10.**(0.03*(idx1-1)+0.02)   ! 1.0e0 to 1.0e3 
      end do

      open(100, file='../Data/Ice_master_fzj_all.dat') 
   
      irec = 0
      do while(.true.) 
        irec = irec + 1 
        if(irec.le.20)  then
          read(100, *)   
        else 
          read(100, fmt=200, end=101) ctemp
200       format(a50)
          if(scan(ctemp, 'Flight').eq.0) then
            backspace(100)
            read(100, fmt=300, end=101) tk, nice, rei
300         format(2x, f6.2, f13.3, f13.3)
            print *, 'irec=', irec, tk, nice, rei
         
            rhi1 = rei
! Ni: 1.0e-4 to 1.0e3
            if(rhi1.ge.1.0e0.and.rhi1.le.1.0e3) then
               idx1=int((log10(rhi1)+0)/0.03)+1
! Ni: 1.0e-1 to 1.0e3
!              if(rhi1.ge.1.0e-1.and.rhi1.le.1.0e3) then
!               idx1=int((log10(rhi1)+1)/0.04)+1
               if(idx1.gt.100.or.idx1.lt.1) then
                 print *, 'idx1 wrong', idx1, rhi1
                 stop
               end if
               idx1 = min(idx1, 100)
               idx1 = max(idx1, 1)
               iT = int(tk - 180)+1
               dist_rh(1, idx1, iT) = dist_rh(1, idx1, iT) + 1.
            end if
          end if  ! end ctemp
        end if ! end irec .le.20
      end do  ! end do while

101   continue
       
      do iT=1,mT 
        total_ncount(1, iT)=sum(dist_rh(1, 1:100, iT))
        if(total_ncount(1, iT).gt.0.1) then
         temp1=total_ncount(1, iT)
         do idx1=1,  100
           dist_rh(1, idx1, iT) = dist_rh(1, idx1, iT) / temp1
         enddo
        else
         print *, 'May be wrong', iT+180 
         print *, 'dist_rh1', dist_rh(1, :, iT)
         dist_rh(1, :, iT) = -9999.99 
        endif
      enddo  ! end iT
      print *, dist_rh(1, 50, 30)

!     write out the distribution  of Ni 

!      cfile1='_annual_clear_cwc1.0e_10.dat'
!      cfile1='_annual_clear_cld01.dat'
!      cfile1='_annual_allsky.dat'
!      cfile1='NEW_Sept0307_clear_cld01.dat'
!      cfile1='_00Oct_cld1.0e-3.dat'
!      cfile1 = '_01annu_cld1.0e-3_T35_60_INCA.dat'
       cfile1 = '_Krameretal2009_data.dat'

      open( 200, file = '../Data/'//                  &
           'REI_pdf'//trim(cfile1))
      write(200, *)  'pdf of RHi'
      write(200, 205)  'REI', 'PDF', 'accum_pdf'
205   format(3x, 3a15)
      do iT=1, mT 
       write(200, 250) 'Temperature (K) ', 180+iT-0.5
       ni_mean = 0.0
       temp1=0.0
       do idx1=1, 100 
         if(dist_rh(1, idx1, iT).gt.0) then
          temp1=temp1 + dist_rh(1, idx1, iT)
          ni_mean = ni_mean + rh(idx1) * dist_rh(1, idx1, iT)
         end if
         accum_pdf(idx1) = temp1
         write(200, 201) rh(idx1), dist_rh(1, idx1, iT), accum_pdf(idx1)
       end do
! Mean
       if(accum_pdf(100).gt.1.0e-5) then
         ni_mean = ni_mean/accum_pdf(100)
       else
         ni_mean = -9999.
       endif

! 25 percentile
       do idx1=1, 100
        idx_25=1
        if(accum_pdf(idx1).gt.0.25) then
          idx_25=idx1
          exit
        end if
       end do
! 50 percentile
       do idx1=1, 100
        idx_50=1
        if(accum_pdf(idx1).gt.0.50) then
          idx_50=idx1
          exit
        end if
       end do
! 75% percentile
       do idx1=1, 100
        idx_75=1
        if(accum_pdf(idx1).gt.0.75) then
          idx_75=idx1
          exit
        end if
       end do

       write(200, 212) 'Tot sample', total_ncount(1, iT)
       write(200, 212) 'Mean', ni_mean
       write(200, 212) '25 percentile',  rh(idx_25), accum_pdf(idx_25)
       write(200, 212) '50 percentile',  rh(idx_50), accum_pdf(idx_50)
       write(200, 212) '75 perdentile',  rh(idx_75), accum_pdf(idx_75)
      end do  ! end iT
201   format(3x, 3e15.4)
211   format(3x, a15, 2e15.4)
212   format(3x, a15, 2e15.8)
210   format(3x, a6, a6, a13, f10.2, a9, a10)
250   format(3x, a15, f10.5)

      end 

      function esw(T)
      real*4 T
      real*4 esw       ! mb
!------------------------------------------------------------------------
!     the formula of saturaftion pressure over water is from WMO, 2000)
!     Reference: WMO, General meteorological standards and recommended parctices, 
!     Appendix A, WMO Technical Regulations, WMO-No. 49, 2000.
!-----------------------------------------------------------------------
      esw=10.79574*(1-273.16/T)                              &
        -5.02800*log10(T/273.16)                            &
        +1.50475*1.0e-4*(1-10**(-8.2969*(T/273.16-1)))      &
        +0.42873*1.0e-3*(10**(-4.76955*(1-273.16/T))-1.0)   &
        +0.78614
      esw=10**esw
      return
      end

      function esi(T) 
      real*4 T
      real*4 esi     ! mb
!------------------------------------------------------------------------------
!    the formula of saturation pressure over ice is from Goff Grathc equation,
!    Smithsonian Tables, 1984
!    Reference: Smithsonian Met. Tables, 5th ed., pp. 350, 1984
!------------------------------------------------------------------------------
      esi=-9.09718*(273.16/T-1.)                          &
         -3.56654*log10(273.16/T)                        &
         +0.876793*(1-T/273.16)                          &
         +log10(6.1071)
      esi=10**esi
      return
      end    

  real*4 function svpl(t)

! using the formular from CAM
! t is the absolute temperature
!
! Water
!
      implicit none
      real*4 t, ts, ps, es
      real*4 e1, e2, f1, f2, f3, f4, f5, f

      ps = 1013.246
      ts = 373.16
      e1 = 11.344*(1.0 - t/ts)
      e2 = -3.49149*(ts/t - 1.0)
      f1 = -7.90298*(ts/t - 1.0)
      f2 = 5.02808*log10(ts/t)
      f3 = -1.3816*(10.0**e1 - 1.0)/10000000.0
      f4 = 8.1328*(10.0**e2 - 1.0)/1000.0
      f5 = log10(ps)
      f  = f1 + f2 + f3 + f4 + f5
      es = (10.0**f)*100.0

      svpl = es/100.   ! mbar

    return
  end function svpl

  real*4 function svpi(t)

     implicit none
!
! Ice
!
      real*4 t, t0, term1, term2, term3, es
      real*4 esl   !++wmh
      real*4 svpl

      t0    = 273.15
      term1 = 2.01889049/(t0/t)
      term2 = 3.56654*log(t0/t)
      term3 = 20.947031*(t0/t)
      es    = 575.185606e10*exp(-(term1 + term2 + term3))

!++wmh
! limit svpi not larger than svpl
      esl = svpl(t)
      svpi = min(esl, es/100.)   ! mbar
!--wmh

    return
  end function svpi
